Quantised motion of an atom in a Gaussian-Laguerre beam 
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We quantise the centre of mass motion of a neutral Cs atom in the presence of a classical 
Gaussian-Laguerreio light field in the large detuning limit. This light field possesses orbital angular 
momentum which is transferred to the atom via spontaneous emissions. We use quantum trajectory 
and analytic methods to solve the master equation for the 2d centre of mass motion with recoil near 
the centre of the beam. For appropriate parameters, we observe heating in both the cartesian and 
polar observables within a few orbits of the atom in the beam. The angular momentum, L, shows a 
rapid diffusion which results in (L) reaching a maximum and then decreasing to zero. We compare 
this with analytic results obtained for an atom illuminated by a superposition of Gaussian-Laguerre 
modes which possess no angular momentum, in the limit of no recoil. 

03.75.Be 



I. INTRODUCTION 

Recently there has been a growing interest, both in 
theory and experiment, in the interaction between mat- 
ter and light fields possessing orbital angular momen- 
tum |jl|. The typical light fields used in these works are 
the Gaussian-Laguerre modes (G-L). Experimental work 
has concentrated on two areas; the study of the motion 
of microscopic particles in Gaussian-Laguerre beams [Q , 
and the use of such beams in atom trapping ||]. The 
theoretical work so far has solely concentrated on the 
semi-classical description of the motion of matter in G- 
L beams M. As neutral atom cooling improves, a truly 
quantum description becomes necessary. Indeed, recent 
experiments in atomic guiding using evanescent waves 
within a hollow optical fibre verge on the quantum regime 
M. In this paper we examine the quantised motion of 
a neutral atom in a Gaussian-Laguerreio mode, taking 
into account the spontaneous recoil while in the far de- 
tuned limit. We are primarily interested in the quan- 
tised transfer of the orbital angular momentum to the 
centre-of-mass (COM) of the atom and how this appears 
in the quantal evolution. In the limit where the inci- 
dent light is far detuned from the atom's natural fre- 
quency, we find that the transfer of angular momentum 
is imparted through the dissipation suffered by the atom 
through spontaneous recoil. We derive a master equa- 
tion for the GOM motion in the limit of large detuning 
for motion near the centre of the beam which includes the 
momentum kick from spontaneous recoil. We select the 
beam and atomic parameters which give an atom that 
can be easily cooled to the lowest COM eigenstate of the 
G-L beam and which suffers a large amount of sponta- 
neous recoil. This results in a rapid increase of angular 
momentum over the orbital time of the atom around the 
center of the beam. The atom, because of the dissipation 
caused by the spontaneous recoil, also heats up. This is 



indicated by an increase in the variances of x and p. 

For the case of a neutral Cesium atom, we solve the 
master equation using quantum trajectories and numeri- 
cally compute expectation values and variances for the 
atom's orbital angular momentum, cartesian position 
and momentum. We solve for an atom initially in a mini- 
mum uncertainty state (1) at the centre of the beam and, 
(2) offset from the centre. We find that the variances in- 
crease with time as the system heats due to the diffusion 
caused by the spontaneous recoil. Surprisingly, we also 
find both in (1) and (2), that the diffusion in the angular 
momentum L, is so rapid that (L) reaches a maximum 
and then decreases to zero. This is to be contrasted with 
the semiclassical work ^ , which shows a general increase 
in the atom's semiclassical angular momentum while in 
the harmonic regime of the light field potential. We fi- 
nally compare this system to a similar master equation 
with no cross terms and no momenta transfer on recoil. 
This corresponds to an atom coupled to a superposition 
of G-L+io and G-L_io modes in the semiclassical limit. 
Here no angular momentum is transferred but the vari- 
ances in both X and y remain much smaller than in the 
case where recoil is included. 



II. GAUSSIAN-LAGUERRE MODES AND 
MASTER EQUATION 



Gaussian Laguerre modes can be produced through 
cylindrical lens mode converters Ui, spiral phase plates 
m and through computer generated holograms M. The 
general form for the electric field for the Gauss-Laguerrejo 
mode is M 
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where w is the beam waist, R is the wavefront radius 
of curvature and $ is the Guoy phase shift. For most 
appUcations the paraxial approximation is quite good. 
In this approximation it has been shown |9|,|l0[| that the 
Guoy phase and Rayleigh curvature are very small for 
Gaussian-Laguerre modes near the beam waist and we 
shall ignore them. We thus restrict the atom's motion to 
a two dimensional plane transverse to the beam propa- 
gation direction, positioned at the beam waist. We shall 
concentrate on the I = +1 mode. Near the centre of the 
beam the exponential dependence on r is small and we 
shall also ignore it. 

When the system is near resonance the internal elec- 
tron degrees of freedom are strongly coupled to the ex- 
ternal spatial degrees of freedom of the atom. To solve 
this fully quantum mechanically, for two spatial degrees 
of freedom, is virtauUy impossible. To simplify matters, 
we work in the far detuned regime. In this regime, the 
master equation, in the dipole approximation, and adia- 
batically eliminating the dynamics of the internal excited 
state (assuming a two- level atom), takes the form |11 
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where Cl is the spatially dependent Rabi frequency, A is 
the detuning, and |z/p = 1 + r2/(2A2). The super op- 
erator Af, describes the effect of a spontaneous emission 
on the atom ||l3| , 
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where $(n) is the dipole emission probability function 
which is defined to be |I3] 



$(n) = — 
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for plane-polarised incident light and 



Hn)^ — {l + {z-nr) , 



(4) 



(5) 



for circularly polarised incident light. Here n is the unit 
vector giving the direction of the spontaneously emitted 
photon, d is the direction of the atomic dipole, and z is 
along the direction of beam propagation. 

Before continuing we will discuss the validity of the 
two assumptions made above. From considerations that 



will follow, we will choose to work with the P- 
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transition in Cs. Firstly, to consider this to be a two-level 
atom we must excite the P+3/2 ^^ '5'+i/2 or P_3/2 ^^ 
S_i/2 transitions. This can be done with incident a± 
circularly polarised light. Plane polarised incident light 
does not result in a closed 2-level atomic system. Thus, 



we must take (^) as the angular probability distribution 
for the direction of spontaneous emission of photons. 



Secondly, we have chosen to work in the dipole ap- 
proximation. We now explain why this approximation is 
still valid even though the COM is quite cold. The elec- 
tronic transitions occur on much faster timescales than 
the dynamics of the COM. Thus, although the COM is 
cooled to where the de-Broglie wavelength of the COM 
degree of freedom is relatively large, the electronic de- 
grees of freedom with which the atom and light interact 
are not cooled and can rapidly follow the slow COM mo- 
tion. Thus the electric dipole approximation holds when 
the spatial gradients of jJlp are small on the scale of the 
electronic clouds involved in the transition and not on 
the scale of the de-Broglie wavelength of the COM. To 
obtain field gradients large enough to significantly ex- 
cite higher multipoles one must be working in the ex- 
treme non-paraxial regime. The degree to which the 
Gauss-Laguerre modes are non-paraxial and thus couple 
to higher multipoles is extremely small. For the phys- 
ical model discussed in this paper we can estimate the 
ratio of the azimuthal spin-orbit force to the transverse 
electric field force to be Fs-i/F^ - 6 x lO^^, Q. The 
omission of higher multipoles is therefore a very good 
approximation for these modes. One consequence of this 
approximation is that the exchange of spin angular mo- 
mentum is conserved separately from the exchange of or- 
bital angular momentum fig ], i.e. there is very little 
spin-orbital coupling. Since $(n) is insensitive to the 
sign of the helicity of the incident light and the Clebsch- 
Gordan coefficients for the transitions in Cs are symmet- 
ric under a -^ —a we are forced to conclude that the 
quantum evolution of the COM is highly insensitive to 
the helicity of the incident light. This conclusion, al- 
though quite reasonable on the atomic scale, seems to be 
at odds with recent experiments on the "classical" mo- 
tion of micro-spheres illuminated by circularly polarised 
Gaussian-Laguerre beams [|]Q|. There are two possible 
physical mechanisms which could generate the observed 
rotation in mesoscopic particles. The transfer of spin an- 
gular momentum to orbital angular momentum could be 
mediated through some bifrefringence present in the ma- 
terial, akin to Beth's original experiment. However, the 
transfer could also be mediated through non-radiative 
off-resonant couplings, i.e. absorption. Although the 
first mechanism might play some small part in the ex- 
periments, recent results seem to indicate that absorp- 
tion is the primary mechanism which couples both the 
spin and orbital angular momentum of the light, to the 
orbital angular momentum of the mesoscopic particles 
p7|,p[. However, in the quantum system studied here, 
the spin-orbit coupling, as we have argued above, for a 
single atom interacting with a Gaussian-Laguerre beam 
should be quite small. An experiment to probe the de- 
pendence of the COM motion of ultra-cold atoms on the 
incident light's helicity would be quite illuminating. 



A. Master Equation 

We begin by recasting the master equation (0) into a 
more convenient form. We first set il = {flo/w) il where 
n — f exp {i9) — X + iy. Next, we obtain the orbital fre- 
quency of the atom in the GLio mode in the harmonic 
approximation to be 
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wavepacket (AP), then the exponential operator in ( [T^ ) 
can be expanded to first order in X. One can then explic- 
itly perform the stochastic average over dn to arrive at 
a Linblad type master equation. This can be converted 
into a c-number Fokker-Planck equation which (since it is 
at most quadratic), will give, in the large energy regime, 
expectation values which agree with the semi-classical 
theory. However, below we will choose the system pa- 
rameters such that /x/3 is large and thus one must solve 
the full quantum dynamics of (O). 



where m is the mass of the atom. We further include the 
following rescalings to give a master equation with co- 
efficients of order unity in the dimensionless quantities, 
X,Y,Ps, and Py, 



X = axX , px = apPx 
y = a^Y , py ^ apPy 
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with T = uigt. This gives [X, Px] = z/3, where /? serves 
as the rescaled Planck's constant. Letting 77 = r/(4A/3), 
we finally obtain 
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For an atom cooled to the recoil limit, PrecoU = ^fc or 
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If the atom is initially in a pure, minimum uncertainty 
state, we have 



APr = ^\ 
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while the effect of the exponential in the superoperator 
Af is to shift the momentum. 
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where fi = kux- We note that the master equation (ll] 
is more involved than the majority of master equations 
previously discussed within the literature. The complica- 
tions here are twofold. In the Linblad form for the master 
equation (pT|), the output channel is not just a momen- 
tum jump (as it is in the near-resonance situation) but 
is a combination of a diffusion in space and momenta. 
Secondly, we note the presence of cross terms between X 
and Y in the dissipation. This cross-coupling frustrates 
the derivation of any simple equation or set of equations 
for the time evolution of ensemble averages. Finally, if 
/i/3 is much smaller than the momentum variance of the 



B. Quantum Trajectories 

The master equation (nTl) must be solved numerically. 
We will use the method of Quantum Trajectories p8[ . 
However, we will be able to solve for the non-unitary 
evolution analytically between jumps. This is possible 
only in the harmonic approximation for motion near the 
centre of the beam. To use this technique we rewrite the 
master equation (hj) in the form 
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(11) and where Cx^y are the x and y components of the recoil 
direction vector n. To apply the method we first choose 
an initial state po = J2i l*i)(^j|- I^ o^'' case we will set 
the initial state to be a pure coherent state. From the 
ensemble po we choose a particular |\l/i) at random and 
apply the following procedure. We generate a uniform 
random variable ^ G [0, 1], and evolve the pure state |^i) 
using the non-unitary Hamiltonian Hnon, 
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i(3 I d'^nC^{n)C{n) = Ho - d{X^ + Y'^) , 
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where 6 = i(3r]. We evolve for a time r^; such that, 
|(*(t^)|*(t(;))P = C- We then apply the "Jump" opera- 
tor C{n), to \^{t(^)), C\^) -^ 1^'), where the vector n is 
generated by two random numbers taken from the distri- 
bution $(z, n), (^). We then renormalise the state |^') 
and repeat the procedure beginning with the generation 
of a new uniform random waiting time, C- At set intervals 
Ts, we store the state. Repeating this whole procedure, 
starting with randomly sampled |^i) from the ensemble 
Po, we generate estimates for the density matrices at the 
times Ts ■ From these estimates we can calculate expecta- 
tions and variances for various observables. In practice, 



unless the final quantities involve the computation of ofF- 
diagonal elements of the p{ts), one only needs to store the 
incremental values for these quantities and not the com- 
plete p{ts) themselves. This produces enormous savings 
on hard memory usage since the storage size of a single 
state l^i), in a truncated Fock state basis of 40, in two 
dimensions, is approximately 25KB in double precision. 
As stated above, in the harmonic approximation, Hnon 
is quadratic and the propagator, C/„o„ = ex.p(iHnonT), 
can be evaluated analytically in the Fock state basis. 
This basis is constructed through the operators 



a^ = (X - iPx)/V2P 



{X + iPx)/^2p , (18) 



for the X and similarly for the y direction. From these 
definitions we can easily see that [aa;,aj,] = 1. In terms 
of these operators the non- unitary Hamiltonian takes the 
form 
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Analytically evaluating the propagator Unon{T) is 
quite tedious and is left to Appendix 1. Inverting, 

l(*»|C^non(T'c)^non('^c)l**)P = C, to obtain r^; is not pos- 
sible analytically and was performed numerically using 
Brent's algorithm [[l9[ . 

To compute the effect of the Jump operator C(n) on 
the state in the Fock basis is quite involved. A num- 
ber of schemes have been used previously to increase 
the efficiency of the quantum trajectory method in the 
case where the quantum jump is a pure momentum kick 
|g^,|l|. These schemes will not work in this case as the 
jump operator, ([iq), is not a pure momentum kick. One 
could perform this jump sequentially by applying the mo- 
mentum shift in the P basis, fast Fourier transforming 
(FFT) into the X basis, and then applying X^. Since 
we are capable of computing the jump analytically we 
doubt that the FFT method would prove more efficient 
in this case. In the more realistic (and more complicated) 
case where we do not make the harmonic approximation, 
no analytical results seem to be possible and the FFT 
method would be necessary. The computation of the ef- 
fect of the Jump operator on the state is given in Ap- 
pendix 2. 

Finally, the generation of the random direction n, for 
the recoil of the atom from the distribution (@) is ex- 
plained in Appendix 3. 



C. Atomic Parameters 

In this section we set out the ingredients necessary for 
a quantum description of the evolution of the COM. We 
require an atom which can be easily cooled to signifi- 
cantly populate the lowest vibrational levels of the G-L 



mode. The dissipation must be large enough to observe 
decay of the system within the timescale of the orbital 
period of the atom around the mode. If the dissipation 
is too weak then the dynamics are washed out and a sec- 
ular approximation, averaging over a cycle time, would 
be a good description of the evolution. This condition 
can be roughly approximated by setting the recoil fre- 
quency equal to the orbital frequency of the harmonic 
well, ujs = hk'^ /2m. 

These two conditions coupled with that of the large de- 
tuning are quite difficult to achieve. For most atoms, the 
"hole" in the beam is very dark and the dissipation coeffi- 
cient 77, is very small. This makes such a mode quite good 
for trapping the atom in a region of low dissipation [BJ . 
To achieve high dissipation 77, we must make F/A as large 
as possible. To enforce the large detuning limit we must 
also have A > F, A > ilo. AH of these conditions imply 
fio ~ r ~ A. For atoms with very low Doppler tem- 
peratures (Calcium), these conditions are hard to meet 
as usually F <C rJo (for modest beam intensities). In- 
stead we choose to look at Cesium, cooled to recoil. For 
the Cs atom, with m ~ 0-665xl0"^^kg, A — 657nm, 
r/27r = 5 • SMHz, with a Gaussian-Laguerre beam of in- 
tensity / = 4W/m^, and width w = 2 x 10~^m, we get 
a Rabi frequency ilo/27r w SMHz. Choosing A '--^ SF 
or A/27r ^ 16 • 4MHz we obtain an orbital frequency of 
uJs/2tt ^ 774Hz. The ratio Erecou/E ground ~ 2 • 6, indi- 
cates a high population in the ground vibrational state 
of the well. Choosing /3 = • 25, we have ax — ■ 6/im, 
Op = • 6 X 10"27^/g^ ^ ^ 0.0125 and p = 2.310. For 
an atom cooled to recoil in a pure minimum uncertainty 
state (|l|) we get AP^ ~ 0-58, and AXr ~ 0-22. 
Since the momentum shift from the bare recoil kick is, 
P ^ P + • 57, and is thus on the same scale as the 
quantum wavepacket we must use the full quantal mas- 
ter equation to describe the evolution. For simplicity, 
we set the initial wavepacket to be a minimum uncer- 
tainty state with equal variances in both P and X so 
that AX = 0-35. In the rescaled time t the period is 
27r//? = 87r w 25. 



D. Numerical Simulation 

We first simulate the dynamics of the atom initially 
positioned at the centre of the trap. We plot only the 
mean and variance for the cartesian position, momen- 
tum and polar angular momentum of the atom. The 
angular momentum operator is given in the position rep- 
resentation as L = {XPy — PjiY). In the Fock basis this 
becomes L = «/3(ajf a^ — a!r,ay-)/2. The code was run on 
two Sun Hyper-Sparcs and took approximately 100 CPU 
hours to produce ~ 300 trajectories. The results of this 
simulation are shown in Figure 1. From the cylindrical 
symmetry we should have (X) = and (Px) = 0. This 
is well approximated in the data as seen in Fig. la. In 
Fig. lb the variances in both X and Y increase with time 



as the dissipation heats the system. We notice a shght 
taihng off of this heating near r = 80. We cannot be 
entirely sure that this is not due to truncation effects as 
the wavepacket spreads out to regions of large laser in- 
tensity and undergoes frequent recoil. The plots of {X) 
and (Y) still remain close to zero and one might suspect 
that this decrease in heating rate may be a real effect. A 
more surprising result can be seen in Fig. Ic. Here, we 
see an initial increase in (L) and (I/^). The variance is at 
all times larger than the average. The average reaches a 
maximum and then decreases towards zero as r increases 
beyond t = 50. This is not a truncation error and is due 
to the very rapid increase in variance after r = 50. The 
diffusion of L becomes so fast that the average quickly 
drops to zero. 

In the second simulation the atom is initially again in 
a minimum uncertainty state but with an initial position 
X = 1, Y = 0, and initial momentum Px = 0, Py = 1. 
In the absence of dissipation the atom will rotate in 
the X — Y plane about the origin with a period of 
Tp = 27r//3 = Stt ~ 25. This second simulation took 
longer to run since there were many more jumps per tra- 
jectory ( to T = 80 ) than in the first simulation. The 
reason for this is that the initial wavepacket is in a re- 
gion of higher laser intensity and thus undergoes more 
frequent spontaneous emissions. Computing times grew 
too lengthy on the Hyper-Sparcs. It became necessary to 
parallelise the code which was then run on a multipro- 
cessor SGI Power Challenger. The results shown below 
consist of approximately 300 trajectories and correspond 
to 50 CPU hours on the Power Challenger. The results 
of the second simulation are shown in Figure 2. We see 
essentially the same overall behavior as in the first simu- 
lation. Heating is shown by the increase in the variances 
of X and Y (Fig. 2b). The variance in L again increases 
faster than the mean and (L), after gaining a maximum, 
falls to zero again (Fig. 2c). In Figure 3 we have plotted 
the X — Y probability distribution of the estimated p at 
various times in the evolution. To gauge the truncation 
errors we have plotted in Fig. 4a a histogram of the num- 
ber of jumps occurring within a given time interval as a 
function of time. From this and the squarish probability 
profile in the X — Y plane we see that truncation effects 
become significant after r = 40. 

In Appendix D we solve the master equation ^^), with- 
out the exponential recoil kick and coupled cross terms in 
the dissipation. This corresponds to semi-classical evo- 
lution (the momentum shift becoming negligible with re- 
spect to the size of the wavepacket) of the atom illumi- 
nated by a coherent superposition of G-L+io and G-L_io 
modes. We see that in this case (AX) increases without 
bound but remains below the values obtained in the nu- 
merical simulations. This is to be expected as there is no 
recoil here. This result gives us a rough analytical check 
on the numerics of the simulation. 



III. CONCLUSION 



In this work we have studied the quantum dynamics of 
the two dimensional spatial COM motion of an ultra-cold 
Cesium atom cooled to recoil moving near the center of a 
Gaussian-Laguerreio light field in the far detuned regime. 
We found that the transfer of orbital angular momentum 
from the light field to the atom is mediated by sponta- 
neous emissions. In the appropriate limit we find that 
the orbital motion of the atom is not very sensitive to 
the helicity of the illuminating light field. We argue that 
the helicity dependence seen in mesoscopic experiments 
may be due to residual birefringence present in the mate- 
rial. For motion near the center of the beam we made an 
harmonic approximation to the light potential. In this 
approximation we solved the quantum dynamical master 
equation using the quantum trajectories method. With 
the harmonic approximation we were able to solve for 
the non-unitary and jump evolutions analytically. This 
significantly increased the computational efficiency of the 
numerics. Two seperate initial conditions were used. In 
the semiclassical regime, in the harmonic approximation 
and with the same initial conditions, the atom gains a 
net angular momentum. In the quantum regime, for the 
particular physical parameters chosen, we found that the 
diffusion of the angular momentum was very high with 
the result that (L) reached a maximum at a time Tc, 
and then decreased to zero. As an analytical check we 
solved for the quantum evolution of the atom illuminated 
by a superposition of G-L+io and G-L_io light in the 
limit of no recoil. The variances found were smaller than 
those seen in the fully quantum computation as expected. 
From similiar studies in other systems one might expect 
that Tc will increase with an increase in the initial state 
energy or a decrease in the dissipation rate. However, 
to explore the parameter space of this system more fully 
and include the complete exponential character of the 
light field potential would require a very significant in- 
crease in computational resources. With the inclusion 
of the exponential in the potential and dissipation terms 
of the master equation we would expect to see shear in 
the evolution of the wavepacket. Partial revivals and cat 
states might be expected as these have been seen in other 
nonlinear potentials. We note that unless some new an- 
alytical techniques are found further investigations into 
this system will be very computationally expensive. 

Finally, experimental probing of (i), for ultra-cold 
atoms in the center of a Gaussian-Laguerre trap could 
be attempted using the azimuthal Doppler shift provided 
by the rotation of the atoms about the axis of the beam 
|25| . One might double the shift by retro-reflecting the 
laser probe beam back through the rotating atomic cloud 
while maintaining the reflected beam's impact parameter 
to the axis to be the same as the incident beam's impact 
parameter. In order to achieve a good Doppler signal it 
may be necessary to work in a regime where the orbital 
frequency of the atomic cloud near the center of the G-L 



mode is higher than that used in this work. 
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we p erform the differentiation on the left hand side of 
( jAq ) and use BCH disentangling to obtain ordinary dif- 
ferential equations relating the gi{e) to the a^. The solu- 
tion of these ODEs for e — I and initial gi (e = 0) gives 
the required normal ordering. 
Differentiating Z-/(e) gives 

^ = nu 

de 

= {g+K+ + .go e^+^'^^+ifo + g- 6^+^^^+ e3^""^^"K^) U 
= [a+K+ + a^K^ + aoKo] U (A7) 



APPENDIX A: MATRIX ELEMENTS FOR Unon 

In this appendix we calculate the matrix elements of 
the non-unitary propagator Unon in the number state ba- 
sis. Although the following calculations are somewhat 
involved the resulting analytic expression proved more 
efficient than a direct numerical integration of the equa- 
tions. We begin with the general Fock state basis matrix 
element 



{m^,my\e "■""""'' /'^\n^,ny) 



e-'^(i-'5)(m^|exp 



{l-5)ala^~ -{a^^ +al) 



X {my\ exp < —IT 



{I - 5)alay - -{af + al) 



(Al) 



where 5 — i(3ri. 

To evaluate this we use Baker-Campbell-Hausdorff 
(BCH) disentangling. The particular methods is due 
to Wilkox |2^. For simplicity we take first the x ma- 
trix elements and drop the x subscript. We then write 
Unon = exp{H) where H = xKq + y{K+ + K-), and 
where Kq, K± are the generators of the Lie group SU(1,1) 
which satisfy 



We now use the commutation properties of SU(1,1) 
and the BCH formula 



e«^ B e-«^ ^B + ^[A,B] + ^ [A, [A, B]] 



^^[A,[A,[A,B]]] 



(A8) 



{A and B are operators while ^ is a scalar) to construct 
the adjoint action table (Table Al) 

Using this information we can evaluate the adjoint ac- 
tions in equation (A7) giving 



{K+ {g+ - 2g+go + e-^s°glg^) + K^ {go - g+e-^'°g-) 
K^ (e"2s".g_)} 



Equating the coefficients of K± , Kq, and re- arranging 
gives 



g- ^a-e^so ^ 

go = flo + g+a- , 

9+ = 5+ (g+a- + 2ao) + a+ 



(A9) 
(AlO) 
(All) 



Ko = ^{a^ a + 1/2) 



K. 



1 



t2 



K- . \a' 



(A2) 



Solving these with the initial condition ^^(e = 0) = and 
setting e = 1 in the final result we get 



[K+,K^]=-Ko , [Ko,K±]^t2Ko 



(A3) 



a+ tan 7 
7 — flQ tan 7 



9- = 



a_ tan 7 
7 — ap tan 7 



(A12) 



-2iT(l — (5) , y ~ +iTS 



(A4) 



To evaluate the matrix element between number states 
we must first re- write Unon in normal order. To do this 
we use the parameter differential methods of Wilkox [g2| . 
Set U{e) = exp [He), this gives 



rf7 



nu 



(A5) 



90 



In 



a-o ■ 
cos 7 sm 7 

7 



(A13) 



where 7^ = a+a- - Oq = (4 + 3(5^ - 8S)t'^. To obtain 
the matrix element of the normally ordered operator, (ie. 
{m\U\n))) we can proceed in a number of ways, however 
the most straightforward is to insert resolutions of unity 
in the coherent state representation, 



and put 7i = a^K^ + a{)K{) + a-K-, [x = oq and 
y = a^ = ajf-). Adopting the ansatz 



{ra\U\n) 



(Pa [(P(3 



■K IT 



'm\a){a\U\l3){(i\n) . (A14) 



U = 6^+^^-'-'^+ e^"^'^-''^'' e^"^'^^ 



K- 



(A6) From (A6) we get 



{m\L(\n) — {a\e 



9+K+ pQoKo pQ-K- 



n 



1 f d^ad^P 



exp \ -gn + i^g+a 



7;9-P' 



offo/2 



-,m 



D 



a*P 



(A15) 



TT^ y \Jn\m\ 



4(|ap + |/3n+eWV/3 



+ 2-9+"*^ + 2-9-'^^ + "^1" + '^2/3* + ^50 



1 /.,2 

exp 



Vnlm! 



(|5- + |5++7i72e^«/^)e^«/4 



(A17) 



where we have used exp(goa^a/2) — 'Y^ (exp (go/2) — 

l)'at'a'Z!-^ The normally ordered element can now be To perform the differentiations in ( |A16| ) we make use of 



expressed as 



the identities 



Ti" J Vn\m\ 



d" 



H..= i-ire^^e 



n!e-^a;"L"(x) = -— (e-=^x"+") , 



dx^ 



d" d™ 



n 



(A16) 



(A18) 



71=0,72=0 



where 



where if„, L" are the Hermite and associated Laguerre 
polynomials. Writing 11 = exp(go/4) exp(aa;^ + by^ + 
cxy)l \Jn\m\ where a — g+/2 ,b = (?-/2 and c = 
exp {go/2) one can, after expanding products, differen- 
tiating and setting 71 = 72 = 0, obtain 



{ra\U\n) 



g9o/4 
y/nlml 



[»/2] 

E 



^0 + fc)aJ;,fc^«-2j(^ - 2fc)!i?2j(0)il2fc(0) 



j^ — mm(0,A) 



2j / I 2fc 



(19) 



where k=j~A, A — {n — rn)/2, and [n/2] represents the integer part of n/2. To construct the matrix elements 
for the combined x ~ y system we have |^) = ^ A{nxTny)\nxTny). Under the non-unitary evolution this becomes 
Unon{T)\^) = I*') = Y. B{rnx,my)\mx,my) where 



B{mx,my) = {■mx,my\Uno7i{T)\ ^ A{nx,ny)\nx, 

Wx ,ny 

= ^ A{nx,ny){vix\U{T)\nx){my\U{T)\ 



yl ' 



(20) 



Denoting the matrix element oiU by Uab, we have the simple relation Bm^,ray — (U ■ A- Vl'^)m^,my- This completes 
the analytic description of the non-unitary propagator. We note that the main numerical overhead occurs in the 
computation of the sum in (n9|) for a given value of r. Lookup tables for the factorial and Hermite functions were 
used to increase efficiency. 



APPENDIX II: MATRIX ELEMENTS FOR THE JUMP OPERATION 



In this appendix we calculate the matrix elements for the Jump operation. To jump we essentially apply the 
operator 



C=ix + iy)e^''('=-*+^''«) , 



(1) 



to the state Ivf) where ex,y are the direction cosines for the recoil kick and fi is related to the momentum transferred. 
We are thus reduced to finding the matrix element 

{mx,my\{x + iy)e'^'^'^^+'yy^\nx,ny) . (2) 

We begin by evaluating {m\ exp{ibx)\n) and then differentiating with respect to b. From section (II) we have x = 
y\572(a^ + a) and putting k = by^P/2 and using the definition of \n) with BCH disentangling we obtain 



(m\e'''^\n) 



1 
\Jn\ra\ 



^mg»KatgiKa^t"|0)e-«'/2 



(3) 



We can pull the annihilation operators to the right using 



m iK.a' __ ina' —iKa' m iKa 



* e*'"°\a + zk)"^ 



Similarly, eyop{iKa)a*'^ — [a* + m)" exp(iKa). We thus must evaluate 



{m\e'^^\n 



-kV2 



VTjJm! 



(a + iK)"(aT + iK)"e^«°|0) 



(4) 



(5) 



Binomially expanding the terms in (0) we finally obtain 






j=o 



n< 



[in) 



n>+n<-2j 



J / ("> -j)! 






n>,-n<],0,-— ) , 



g{m,n,b) , 



(6) 

(7) 



where n> = max(n,77i),n< = min(n, tti) and T is the generalised hypergeometric function |23]. Now to obtain the 
matrix element (m|e*''^|n) we differentiate with respect to k to get 



__ Ijl /n>!e-«/2(m) 



-K*F([-n>,-n<],0, 5-) 



2 y n<! K3n>! 

+ Ki(n + m)F([-n>,-n<],0, ^) + 2n>n<F([-n> + l,-n< + 1],0, 

= J-'{m, n, 6) . 

The complete two dimensional matrix element may now be expressed as 

{nix ,my\{x + iy)e'i'^'''^+'yy^ \n^, Uy) 

= T{mx,nj:, fiex)G{my,ny, ^ey) + iG{mx, rix, fJ-£x)J^imy, ny,fiey) 
= T{mx,my]nx,ny) . 



(8) 
(9) 



(10) 



Denoting the pre-jump state by |^) = ^ A{nxTny)\nx,ny), then |(\I'|\I')|^ < 1 after then jump. Letting the state 
after the jump be \^)j = C\^) = ^ B{mx,my)\mx,my) then the coefficients B{mx,my) are given by 



B{mx,my) ^ ^ A{nx,ny)T{mx,my;nx,ny) 



(11) 



The main numerical overhead occurs in the calculation of Q and J-'. Since the jump directions are chosen randomly, 
the values for these functions must be computed each time. However, lookup tables for the generalised hypergeometric 
function are used and only two evaluations of F are needed for each jump computation. After a jump the state is 
renormalised. 



APPENDIX III: GENERATION OF RANDOM 
EMISSION VECTORS 

To generate the random direction vectors, n{(j), 9), sam- 
pled from the probability distribution (0) we first note 
that the probability is independent of the angle (j). Thus 
we choose to be a uniform random variable (j) € [0, 27r]. 
To generate the 6, a random variable sampled from the 



probability distribution 1 + cos^ 9, we use the method of 
cumulative inversion. We set the cumulative probability 
in 9 equal to a uniform random variable and invert to 
obtain 9, 



ain9d9 (1 + cos^ 61) 

IGtt 



(1) 



gives 



where «, = 2 — 4e + VS — 16e + 16e^. 



(2) 



APPENDIX IV: CALDEIRA-LEGGETT MODEL 

The dynamics of an atom moving in the semiclassi- 
cal regime (/i = 0), near the beam axis for a beam in a 
superposition of the G-L+10 and G-L_io is that of the 
Caldeira-Leggett model J2^]. This mode does not pos- 
sess any orbital angular momentum and thus there are 
no cross terms in the dissipation in the quantum master 
equation. The equation separates into x and y compo- 
nents. In the position representation the x-component 
equation is 



where Xq = X -\- X' and C,q = X — X' , the coordinates 
for the initial state. Transforming to the barred variables 
we have a = af]. Collecting all of the above we get 

m 



a^ cos^ V -t- ?7/3^ {v - 7^ sin 2v\ -f -^^ sin^ v , (9) 



4ct^ 



where v = tout = (3t = t. Plotting this in Figure 6 for the 
initial conditions fj = 0.0125, (3 = .25 and a = 0.35 we 
see that the variance generally increases with time r. It 
remains smaller than the variance calculated numerically 
for the G-Lio mode with recoil. 



<- 



- V (« - aJ) + li.^'- *") - "iCii - ^") 



In Gm , the master equation, in the limits 
7^0 , r — > 00 , ks^T -^ ry 



IS 



'jtP = 



-\{^l^^l,) + ^{x^-x'^)^^r^{x-x'f 



(1) 
(2) 

P ■ 
(3) 



To convert between these two equations we set r = /3i, 
X — (3x and 77 = f](3^. This gives uj\ = ffi and variance, 
(^2) = /32(j;2). In the limit (|), the coefficients (in [|4), 
equations 4.41, 4.42, 4.43) become 

'7 



[v-sa\2v] , (4) 

[sin V — V cos lA , (5) 

(6) 

„ ^ L = N = Tr^4^^ and where 

I cot V ' 2 sin V J — I 

V = ujRt. We begin with the density matrix (6.14) in |2J] 
we can use their solution for the evolved p{t) in (6.15), 



A 



B 



C 



2lur sin v 

^ 

uiR sin" ly ^ 

V 
2a; 7? sin 1/ 



V sin 2v 

2 



where K 



K 



t^R 



p[t) = iV2 exp {-ACl - BXl + iQiXiC) 



(7) 



where Xi = X + X' , (i — X ~ X' , are the sum and 
relative coordinates for the evolved density matrix. Now 
with the initial condition (P) = , we have (cc^) — gg 
for any t. So we just proceed to carefully evaluate the co- 
efficient B. Beginning with an initial state in the barred 
variables 



PA = 



V2: 



na^ 



(8) 
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TABLE L Table of adjoint action of SU(1,1): e""^^ Ffe e" 





r^rfc 


Ka 


K+ 


K^ 


Ka 


Ka 


e'-K^ 


e-^-K_ 


K+ 


Ka - 2xK+ 


K+ 


K- - xKo + x''K^ 


K^ 


Ko + 2xK. 


K+ +xKa + x^K- 


K^ 
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FIG. 1. Plots of (a) (X), (Y), (b) (AX), (AY), (X solid, 
Y dashed), and (c) (L), (AL), (solid, dashed), for Simu- 
lation 1 where the initial state is a minimum uncertainty 
state at the origin. 




_ FIG. 2. Plots of (a) jX), (Y), (b) (AX), {AY}, (X solid, 
Y dashed), and (c) (L), (AL), (solid, dashed), for Simu- 
lation 2 where the initial state is a minimum uncertainty 
state offset from the origin. 
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FIG. 3. Plots (a) to (g) are the X — Y probability densities for the estimated density matrix in Simulation 2 at the times r = 4mr 
where n — 0, ... 6. 
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(b) 




20 40 X 60 80 100 



FIG. 4. Plot (a) is the histograni of jump frequencies for Simulation 2, i.e. numer of jumps which occurred 
between time t and r + 5, while (b) plots (AX(t)) as an analytic comparison to Simulation 1. Here we have 
neglected the cross terms and recoil momentum kick in the master equation. This is solved analytically. See 
Appendix D for details. 
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